{r} [Zrodlo] danych] (http://www.cs.put.poznan.pl/dbrzezinski/teaching/zed/wuhan_blood_sample_data_Jan_Feb_2020.xlsx) {r} [Opis danych] (https://www.nature.com/articles/s42256-020-0180-7)
raport powinien zaczynać się od rozdziału podsumowującego całą analizę, streszczającego najważniejsze spostrzeżenia analityka
Kod wyliczający wykorzystane biblioteki. Kod zapewniający powtarzalność wyników przy każdym uruchomieniu raportu na tych samych danych. Kod pozwalający wczytać dane z pliku. Kod czyszczący dane (np. zmiany nazw kolumn, tranformacja danych do innych jednostek, decyzje dotyczące brakujących danych). Sekcję podsumowującą rozmiar zbioru i podstawowe statystyki. Szczegółową analizę wartości atrybutów. Sekcję sprawdzającą korelacje między zmiennymi; sekcja ta powinna zawierać jakąś formę graficznej prezentacji korelacji. Interaktywny wykres lub animację prezentującą zmianę wybranych atrybutów w czasie. Sekcję próbującą stworzyć klasyfikator przewidujący czy dany pacjent przeżyje (w tej sekcji należy wykorzystać wiedzę z pozostałych punktów oraz wykonać dodatkowe czynności, które mogą poprawić trafność predykcji); dobór parametrów modelu oraz oszacowanie jego skuteczności powinny zostać wykonane za pomocą techniki podziału zbioru na dane uczące, walidujące i testowe; trafność klasyfikacji powinna zostać oszacowana na podstawie kliku wybranych (i uzasadnionych) miar oceny klasyfikacji. Analizę ważności atrybutów najlepszego znalezionego modelu. Dodatkowo punktowane będzie wykonanie analizy typowej dla danych klinicznych, np. regresji logistycznej wraz z wzięciem pod uwagę czynników zakłócających (ang. confounding factors) lub regresji Coxa (ang. Cox Proportional-Hazards Model).
Analityk nie musi, a nawet nie powinien, ograniczać się do powyższych punktów. Wszelkie dodatkowe techniki analizy danych, wizualizacje, spostrzeżenia będą pozytywnie wpływały na ocenę.
Ewentualne konkluzje, znalezione zależności warto potwierdzić dokonując sprawdzenia istniejących wyników badań w literaturze naukowej (np. na Google Scholar czy PubMed).
TBA
library(xlsx)
library(DT)
library(knitr)
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(tidyr)
library(janitor)##
## Attaching package: 'janitor'## The following objects are masked from 'package:stats':
##
## chisq.test, fisher.testlibrary(imputeTS)## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoolibrary(lares)##
## Attaching package: 'lares'## The following object is masked from 'package:janitor':
##
## crosstablibrary(plotly)## Loading required package: ggplot2##
## Attaching package: 'plotly'## The following object is masked from 'package:ggplot2':
##
## last_plot## The following object is masked from 'package:stats':
##
## filter## The following object is masked from 'package:graphics':
##
## layoutfilename <- 'res/wuhan_blood_sample_data_Jan_Feb_2020.xlsx'
raw_data <- read.xlsx(filename, 1)
raw_data <- as_tibble(raw_data)# OK
#how many NA in variables
mean(is.na(raw_data %>% select(Hypersensitive.cardiac.troponinI : RBC.distribution.width.SD)))## [1] 0.8770596# var min i max
min(raw_data %>% select(Hypersensitive.cardiac.troponinI : RBC.distribution.width.SD), na.rm = TRUE)## [1] -1max(raw_data %>% select(Hypersensitive.cardiac.troponinI : RBC.distribution.width.SD), na.rm = TRUE)## [1] 50000glimpse(raw_data)## Rows: 6,120
## Columns: 81
## $ PATIENT_ID <dbl> 1, NA...
## $ RE_DATE <dttm> 2020...
## $ age <dbl> 73, 7...
## $ gender <dbl> 1, 1,...
## $ Admission.time <dttm> 2020...
## $ Discharge.time <dttm> 2020...
## $ outcome <dbl> 0, 0,...
## $ Hypersensitive.cardiac.troponinI <dbl> NA, N...
## $ hemoglobin <dbl> NA, 1...
## $ Serum.chloride <dbl> NA, N...
## $ Prothrombin.time <dbl> NA, N...
## $ procalcitonin <dbl> NA, N...
## $ eosinophils... <dbl> NA, 0...
## $ Interleukin.2.receptor <dbl> NA, N...
## $ Alkaline.phosphatase <dbl> NA, N...
## $ albumin <dbl> NA, N...
## $ basophil... <dbl> NA, 0...
## $ Interleukin.10 <dbl> NA, N...
## $ Total.bilirubin <dbl> NA, N...
## $ Platelet.count <dbl> NA, 1...
## $ monocytes... <dbl> NA, 1...
## $ antithrombin <dbl> NA, N...
## $ Interleukin.8 <dbl> NA, N...
## $ indirect.bilirubin <dbl> NA, N...
## $ Red.blood.cell.distribution.width. <dbl> NA, 1...
## $ neutrophils... <dbl> NA, 6...
## $ total.protein <dbl> NA, N...
## $ Quantification.of.Treponema.pallidum.antibodies <dbl> NA, N...
## $ Prothrombin.activity <dbl> NA, N...
## $ HBsAg <dbl> NA, N...
## $ mean.corpuscular.volume <dbl> NA, 9...
## $ hematocrit <dbl> NA, 3...
## $ White.blood.cell.count <dbl> NA, 3...
## $ Tumor.necrosis.factorÎ. <dbl> NA, N...
## $ mean.corpuscular.hemoglobin.concentration <dbl> NA, 3...
## $ fibrinogen <dbl> NA, N...
## $ Interleukin.1β <dbl> NA, N...
## $ Urea <dbl> NA, N...
## $ lymphocyte.count <dbl> NA, 0...
## $ PH.value <dbl> 7.415...
## $ Red.blood.cell.count <dbl> NA, 4...
## $ Eosinophil.count <dbl> NA, 0...
## $ Corrected.calcium <dbl> NA, N...
## $ Serum.potassium <dbl> NA, N...
## $ glucose <dbl> NA, N...
## $ neutrophils.count <dbl> NA, 2...
## $ Direct.bilirubin <dbl> NA, N...
## $ Mean.platelet.volume <dbl> NA, 1...
## $ ferritin <dbl> NA, N...
## $ RBC.distribution.width.SD <dbl> NA, 4...
## $ Thrombin.time <dbl> NA, N...
## $ X...lymphocyte <dbl> NA, 2...
## $ HCV.antibody.quantification <dbl> NA, N...
## $ D.D.dimer <dbl> NA, N...
## $ Total.cholesterol <dbl> NA, N...
## $ aspartate.aminotransferase <dbl> NA, N...
## $ Uric.acid <dbl> NA, N...
## $ HCO3. <dbl> NA, N...
## $ calcium <dbl> NA, N...
## $ Amino.terminal.brain.natriuretic.peptide.precursor.NT.proBNP. <dbl> NA, N...
## $ Lactate.dehydrogenase <dbl> NA, N...
## $ platelet.large.cell.ratio. <dbl> NA, 3...
## $ Interleukin.6 <dbl> NA, N...
## $ Fibrin.degradation.products <dbl> NA, N...
## $ monocytes.count <dbl> NA, 0...
## $ PLT.distribution.width <dbl> NA, 1...
## $ globulin <dbl> NA, N...
## $ γ.glutamyl.transpeptidase <dbl> NA, N...
## $ International.standard.ratio <dbl> NA, N...
## $ basophil.count... <dbl> NA, 0...
## $ X2019.nCoV.nucleic.acid.detection <dbl> NA, N...
## $ mean.corpuscular.hemoglobin. <dbl> NA, 3...
## $ Activation.of.partial.thromboplastin.time <dbl> NA, N...
## $ High.sensitivity.C.reactive.protein <dbl> NA, N...
## $ HIV.antibody.quantification <dbl> NA, N...
## $ serum.sodium <dbl> NA, N...
## $ thrombocytocrit <dbl> NA, 0...
## $ ESR <dbl> NA, N...
## $ glutamic.pyruvic.transaminase <dbl> NA, N...
## $ eGFR <dbl> NA, N...
## $ creatinine <dbl> NA, N...mean(is.na(raw_data %>% select(Hypersensitive.cardiac.troponinI : RBC.distribution.width.SD)))## [1] 0.8770596(np. zmiany nazw kolumn, tranformacja danych do innych jednostek, decyzje dotyczące brakujących danych)
# NOT OK
#ogarnac te brakujace daty
#replace -1 with NA
raw_data[raw_data==-1]<-NA
#PATIENT_ID
id_filled <- raw_data %>% fill(PATIENT_ID)
#OTHER NA VALUES
#is there any rows that doesn't have NA value?
full_row_count <-dim(id_filled[complete.cases(id_filled),])[1]
print(full_row_count)## [1] 0# NA values in variables:
mean(is.na(id_filled %>% select(Hypersensitive.cardiac.troponinI : RBC.distribution.width.SD)))## [1] 0.8770748# there are two patients with negative days, as their only blood sample results arrived one day after their clinical outcome
#remove rows where all variables are empty
vars <- colnames(id_filled)[-(1:7)]
no_empty_rows<- id_filled[rowSums(is.na(id_filled[vars])) != length(vars), ]
# NA values in variables:
mean(is.na(no_empty_rows %>% select(Hypersensitive.cardiac.troponinI : RBC.distribution.width.SD)))## [1] 0.8658041colnames_cleaned <- no_empty_rows %>% clean_names()NA Values:
# NOT OK
## ten kod daje oczekiwany wynik, ale z pewnoscia da sie (trzeba) go przerobic
### czy nazwa R wziela sie od slowa "rozpacz"?
clean_NA_in_column<-function(column){
not_NA_count<-sum(!is.na(column))
if (not_NA_count>=2){
column <- na_interpolation(column, option = "linear")
column
}
else if (not_NA_count==1){
val <- first(na.omit(column))
column[is.na(column)] <- val
column
}
column
}
clean_NA <- function(input){
res <- sapply(input, function(x) clean_NA_in_column(x))
res
}
vars <- colnames(colnames_cleaned)[-(1:7)]
colnames_cleaned_dim<- dim(colnames_cleaned)
cleaned <- data.frame(matrix(ncol=81,nrow=0))
for (p_id in distinct(colnames_cleaned, patient_id)$patient_id){
#for (p_id in 1:10){
patient_data <- colnames_cleaned%>%filter(patient_id==p_id)
row_count<- dim(patient_data)[1]
dt=data.frame(matrix(ncol=0,nrow=row_count))
for (var in vars){
column <- patient_data%>%select(var)
res<-clean_NA_in_column(column)
dt<-cbind(dt,res)
}
patient_data[-(1:7)] = dt
cleaned<-rbind(cleaned, patient_data)
}## Note: Using an external vector in selections is ambiguous.
## i Use `all_of(var)` instead of `var` to silence this message.
## i See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
## This message is displayed once per session.“Sekcja podsumowująca rozmiar zbioru i podstawowe statystyki” Szczegółowa analiza wartości atrybutów.
DT::datatable(cleaned, style="bootstrap", filter = "top", rownames = FALSE, extensions = "Buttons", options = list(scrollX = TRUE, dom = 'Bfrtip', buttons = c('copy', 'csv', 'excel')))## Warning in instance$preRenderHook(instance): It seems your data is too big
## for client-side DataTables. You may consider server-side processing: https://
## rstudio.github.io/DT/server.htmlall_dim <-dim(cleaned)
tmp <- cleaned %>% select(patient_id, outcome, gender) %>% group_by(patient_id) %>% summarise(outcome_count = mean(outcome), gender_count = mean(gender))## `summarise()` ungrouping output (override with `.groups` argument)# ile pacjentow
patients_count <- length(distinct(cleaned, patient_id)$patient_id)
# ile pomiarow
measurements_count <- all_dim[1]
# ile srednio pomiarow na pacjenta
mean_measures <- measurements_count/patients_count
# ile kobiet/mezczyzn
genders <- tmp %>% count(gender_count) #Male 224, Female 151 # WYKRES
# ile umarlo/przezylo
outcomes <- tmp %>% count(outcome_count) #201 recovered, 174 died # WYKRES
# ile kolumn
columns_count <- all_dim[2]
# ile atrybutow
vars_count <- all_dim[2]-7
# ile zostalo wartosci NA
na_left <- round(100*mean(is.na(cleaned)), 0)
titles <- c("Liczba pacjentów", "Liczba pomiarów", "Średnia liczba pomiarów", "K", "M", "Smierc", "Zycie :(", "Liczba wierszy", "Liczba zmiennych", "Brakujace wartosci [%]")
values <- c(patients_count,measurements_count,mean_measures,1, 2, 0,1, columns_count,vars_count, na_left)
info_table <- data_frame(
titles,
values
)## Warning: `data_frame()` is deprecated as of tibble 1.1.0.
## Please use `tibble()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.knitr::kable(info_table)| titles | values |
|---|---|
| Liczba pacjentów | 361.00000 |
| Liczba pomiarów | 5606.00000 |
| Srednia liczba pomiarów | 15.52909 |
| K | 1.00000 |
| M | 2.00000 |
| Smierc | 0.00000 |
| Zycie :( | 1.00000 |
| Liczba wierszy | 81.00000 |
| Liczba zmiennych | 74.00000 |
| Brakujace wartosci [%] | 8.00000 |
#wykres death_bothgenders, alive_bothgenders
invisible(remove(tmp))
plot_data <- cleaned%>%select(patient_id, gender, admission_time, discharge_time, outcome) %>%group_by(patient_id) %>% summarise_all(funs(first))## Warning: `funs()` is deprecated as of dplyr 0.8.0.
## Please use a list of either functions or lambdas:
##
## # Simple named list:
## list(mean = mean, median = median)
##
## # Auto named with `tibble::lst()`:
## tibble::lst(mean, median)
##
## # Using lambdas
## list(~ mean(., trim = .2), ~ median(., na.rm = TRUE))
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.admission_plot <- ggplot() +
coord_cartesian() +
scale_color_hue() +
facet_wrap(~outcome) +
layer(data=plot_data,
mapping=aes(
x=admission_time,
y=discharge_time, colour=factor(gender)
),
stat="identity",
geom="point",
position=
position_jitter()
)
ggplotly(admission_plot)## Warning: `group_by_()` is deprecated as of dplyr 0.7.0.
## Please use `group_by()` instead.
## See vignette('programming') for more help
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.Szczegółową analizę wartości atrybutów. TBA
“Sekcję sprawdzającą korelacje między zmiennymi; sekcja ta powinna zawierać jakąś formę graficznej prezentacji korelacji.”
# OK
#TODO wiecej wynikow
var_df<-cleaned[-(1:7)]
plot_data <- corr_cross(var_df,
max_pvalue = 0.05,
top = 30,
rm.na=TRUE)## Returning only the top 30. You may override with the 'top' argument## Warning in theme_lares2(legend = "top"): Font Arial Narrow is not installed, has
## other name, or can't be found#plot(plot_data)
ggplotly(plot_data)Interaktywny wykres lub animację prezentującą zmianę wybranych atrybutów w czasie.
timeline_plot <- ggplot() +
coord_cartesian() +
scale_color_hue() +
layer(data=cleaned,
mapping=aes(
x=re_date,
y=hypersensitive_cardiac_troponin_i, group=patient_id
),
stat="identity",
geom="point",
position=
position_jitter()
)
ggplotly(timeline_plot)timeline_plot <- ggplot(cleaned, aes(x=re_date, y=serum_chloride, colour=factor(patient_id), group=patient_id)) + geom_line() + geom_point() + facet_wrap(~outcome)
ggplotly(timeline_plot)#timeline_data <-cleaned%>%select(patient_id,re_date, hemoglobin, glucose) %>% transmute(id=patient_id, re_date=as.Date(re_date), hemoglobin, glucose) %>%group_by(id, re_date) %>% summarise(avg_hemoglobin = mean(hemoglobin), avg_glucose=mean(glucose)) %>%group_by(id) %>% mutate(day_no=re_date-min(re_date)+1, avg_hemoglobin=round(avg_hemoglobin,2), avg_glucose=round(avg_glucose,2))%>%ungroup() %>%select(day_no,avg_hemoglobin, avg_glucose)
timeline_data <-cleaned%>%select(patient_id,re_date, hemoglobin, glucose, outcome, gender) %>% transmute(id=patient_id, re_date=as.Date(re_date), hemoglobin, glucose, outcome, gender) %>%group_by(id, re_date) %>% summarise(avg_hemoglobin = mean(hemoglobin), avg_glucose=mean(glucose), outcome=first(outcome), gender=first(gender)) %>%group_by(id) %>% mutate(day_no=re_date-min(re_date)+1, avg_hemoglobin=round(avg_hemoglobin,2), avg_glucose=round(avg_glucose,2))%>%ungroup() %>%select(day_no,avg_hemoglobin, avg_glucose, outcome, gender)## `summarise()` regrouping output by 'id' (override with `.groups` argument)timeline_plot <- ggplot(timeline_data,aes(day_no,avg_hemoglobin, shape=factor(gender)))+geom_point(color="blue")+geom_point(aes(day_no,avg_glucose, shape=factor(gender)),color="black") + facet_wrap(~outcome)
ggplotly(timeline_plot)## Don't know how to automatically pick scale for object of type difftime. Defaulting to continuous.klasyfikator przewidujący czy dany pacjent przeżyje (w tej sekcji należy wykorzystać wiedzę z pozostałych punktów oraz wykonać dodatkowe czynności, które mogą poprawić trafność predykcji); dobór parametrów modelu oraz oszacowanie jego skuteczności powinny zostać wykonane za pomocą techniki podziału zbioru na dane uczące, walidujące i testowe; trafność klasyfikacji powinna zostać oszacowana na podstawie kliku wybranych (i uzasadnionych) miar oceny klasyfikacji.
analiza typowa dla danych klinicznych, np.: